{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc13e445",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e9872377",
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'Marketret_mon_stock2023.csv'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[1], line 35\u001b[0m\n\u001b[0;32m     33\u001b[0m pd\u001b[38;5;241m.\u001b[39mset_option(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mdisplay.max_columns\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m)\n\u001b[0;32m     34\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mpandas\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtseries\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01moffsets\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m MonthEnd \u001b[38;5;66;03m# 月末\u001b[39;00m\n\u001b[1;32m---> 35\u001b[0m Market_ret \u001b[38;5;241m=\u001b[39m pd\u001b[38;5;241m.\u001b[39mread_csv(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mMarketret_mon_stock2023.csv\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m     36\u001b[0m Market_ret[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mmonth\u001b[39m\u001b[38;5;124m'\u001b[39m] \u001b[38;5;241m=\u001b[39m pd\u001b[38;5;241m.\u001b[39mto_datetime(Market_ret[\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mmonth\u001b[39m\u001b[38;5;124m'\u001b[39m], \u001b[38;5;28mformat\u001b[39m\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m%\u001b[39m\u001b[38;5;124mb \u001b[39m\u001b[38;5;124m%\u001b[39m\u001b[38;5;124mY\u001b[39m\u001b[38;5;124m'\u001b[39m) \u001b[38;5;241m+\u001b[39m MonthEnd(\u001b[38;5;241m0\u001b[39m)\n\u001b[0;32m     37\u001b[0m Market_ret\u001b[38;5;241m.\u001b[39mset_index(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mmonth\u001b[39m\u001b[38;5;124m'\u001b[39m, inplace\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n",
      "File \u001b[1;32mc:\\Users\\ma970\\anaconda3\\Lib\\site-packages\\pandas\\io\\parsers\\readers.py:1026\u001b[0m, in \u001b[0;36mread_csv\u001b[1;34m(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, date_format, dayfirst, cache_dates, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, encoding_errors, dialect, on_bad_lines, delim_whitespace, low_memory, memory_map, float_precision, storage_options, dtype_backend)\u001b[0m\n\u001b[0;32m   1013\u001b[0m kwds_defaults \u001b[38;5;241m=\u001b[39m _refine_defaults_read(\n\u001b[0;32m   1014\u001b[0m     dialect,\n\u001b[0;32m   1015\u001b[0m     delimiter,\n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m   1022\u001b[0m     dtype_backend\u001b[38;5;241m=\u001b[39mdtype_backend,\n\u001b[0;32m   1023\u001b[0m )\n\u001b[0;32m   1024\u001b[0m kwds\u001b[38;5;241m.\u001b[39mupdate(kwds_defaults)\n\u001b[1;32m-> 1026\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m _read(filepath_or_buffer, kwds)\n",
      "File \u001b[1;32mc:\\Users\\ma970\\anaconda3\\Lib\\site-packages\\pandas\\io\\parsers\\readers.py:620\u001b[0m, in \u001b[0;36m_read\u001b[1;34m(filepath_or_buffer, kwds)\u001b[0m\n\u001b[0;32m    617\u001b[0m _validate_names(kwds\u001b[38;5;241m.\u001b[39mget(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mnames\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m))\n\u001b[0;32m    619\u001b[0m \u001b[38;5;66;03m# Create the parser.\u001b[39;00m\n\u001b[1;32m--> 620\u001b[0m parser \u001b[38;5;241m=\u001b[39m TextFileReader(filepath_or_buffer, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwds)\n\u001b[0;32m    622\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m chunksize \u001b[38;5;129;01mor\u001b[39;00m iterator:\n\u001b[0;32m    623\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m parser\n",
      "File \u001b[1;32mc:\\Users\\ma970\\anaconda3\\Lib\\site-packages\\pandas\\io\\parsers\\readers.py:1620\u001b[0m, in \u001b[0;36mTextFileReader.__init__\u001b[1;34m(self, f, engine, **kwds)\u001b[0m\n\u001b[0;32m   1617\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39moptions[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mhas_index_names\u001b[39m\u001b[38;5;124m\"\u001b[39m] \u001b[38;5;241m=\u001b[39m kwds[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mhas_index_names\u001b[39m\u001b[38;5;124m\"\u001b[39m]\n\u001b[0;32m   1619\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhandles: IOHandles \u001b[38;5;241m|\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[1;32m-> 1620\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_engine \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_make_engine(f, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mengine)\n",
      "File \u001b[1;32mc:\\Users\\ma970\\anaconda3\\Lib\\site-packages\\pandas\\io\\parsers\\readers.py:1880\u001b[0m, in \u001b[0;36mTextFileReader._make_engine\u001b[1;34m(self, f, engine)\u001b[0m\n\u001b[0;32m   1878\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mb\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m mode:\n\u001b[0;32m   1879\u001b[0m         mode \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mb\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m-> 1880\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhandles \u001b[38;5;241m=\u001b[39m get_handle(\n\u001b[0;32m   1881\u001b[0m     f,\n\u001b[0;32m   1882\u001b[0m     mode,\n\u001b[0;32m   1883\u001b[0m     encoding\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39moptions\u001b[38;5;241m.\u001b[39mget(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mencoding\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m),\n\u001b[0;32m   1884\u001b[0m     compression\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39moptions\u001b[38;5;241m.\u001b[39mget(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcompression\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m),\n\u001b[0;32m   1885\u001b[0m     memory_map\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39moptions\u001b[38;5;241m.\u001b[39mget(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmemory_map\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28;01mFalse\u001b[39;00m),\n\u001b[0;32m   1886\u001b[0m     is_text\u001b[38;5;241m=\u001b[39mis_text,\n\u001b[0;32m   1887\u001b[0m     errors\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39moptions\u001b[38;5;241m.\u001b[39mget(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mencoding_errors\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mstrict\u001b[39m\u001b[38;5;124m\"\u001b[39m),\n\u001b[0;32m   1888\u001b[0m     storage_options\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39moptions\u001b[38;5;241m.\u001b[39mget(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mstorage_options\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m),\n\u001b[0;32m   1889\u001b[0m )\n\u001b[0;32m   1890\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhandles \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[0;32m   1891\u001b[0m f \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mhandles\u001b[38;5;241m.\u001b[39mhandle\n",
      "File \u001b[1;32mc:\\Users\\ma970\\anaconda3\\Lib\\site-packages\\pandas\\io\\common.py:873\u001b[0m, in \u001b[0;36mget_handle\u001b[1;34m(path_or_buf, mode, encoding, compression, memory_map, is_text, errors, storage_options)\u001b[0m\n\u001b[0;32m    868\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(handle, \u001b[38;5;28mstr\u001b[39m):\n\u001b[0;32m    869\u001b[0m     \u001b[38;5;66;03m# Check whether the filename is to be opened in binary mode.\u001b[39;00m\n\u001b[0;32m    870\u001b[0m     \u001b[38;5;66;03m# Binary mode does not support 'encoding' and 'newline'.\u001b[39;00m\n\u001b[0;32m    871\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m ioargs\u001b[38;5;241m.\u001b[39mencoding \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mb\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m ioargs\u001b[38;5;241m.\u001b[39mmode:\n\u001b[0;32m    872\u001b[0m         \u001b[38;5;66;03m# Encoding\u001b[39;00m\n\u001b[1;32m--> 873\u001b[0m         handle \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mopen\u001b[39m(\n\u001b[0;32m    874\u001b[0m             handle,\n\u001b[0;32m    875\u001b[0m             ioargs\u001b[38;5;241m.\u001b[39mmode,\n\u001b[0;32m    876\u001b[0m             encoding\u001b[38;5;241m=\u001b[39mioargs\u001b[38;5;241m.\u001b[39mencoding,\n\u001b[0;32m    877\u001b[0m             errors\u001b[38;5;241m=\u001b[39merrors,\n\u001b[0;32m    878\u001b[0m             newline\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[0;32m    879\u001b[0m         )\n\u001b[0;32m    880\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m    881\u001b[0m         \u001b[38;5;66;03m# Binary mode\u001b[39;00m\n\u001b[0;32m    882\u001b[0m         handle \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mopen\u001b[39m(handle, ioargs\u001b[38;5;241m.\u001b[39mmode)\n",
      "\u001b[1;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'Marketret_mon_stock2023.csv'"
     ]
    }
   ],
   "source": [
    "import numpy as np # 数据处理最重要的模块\n",
    "import pandas as pd # 数据处理最重要的模块\n",
    "import scipy.stats as stats # 统计模块\n",
    "import scipy\n",
    "# import pymysql  # 导入数据库模块\n",
    "\n",
    "from datetime import datetime # 时间模块\n",
    "import statsmodels.formula.api as smf  # OLS regression\n",
    "\n",
    "# import pyreadr # read RDS file\n",
    "\n",
    "from matplotlib import style\n",
    "import matplotlib.pyplot as plt  # 画图模块\n",
    "import matplotlib.dates as mdates\n",
    "\n",
    "from matplotlib.font_manager import FontProperties # 作图中文\n",
    "from pylab import mpl\n",
    "#mpl.rcParams['font.sans-serif'] = ['SimHei']\n",
    "#plt.rcParams['font.family'] = 'Times New Roman'\n",
    "\n",
    "#输出矢量图 渲染矢量图\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'svg'\n",
    "\n",
    "from IPython.core.interactiveshell import InteractiveShell # jupyter运行输出的模块\n",
    "#显示每一个运行结果\n",
    "InteractiveShell.ast_node_interactivity = 'all'\n",
    "\n",
    "#设置行不限制数量\n",
    "#pd.set_option('display.max_rows',None)\n",
    "\n",
    "#设置列不限制数量\n",
    "pd.set_option('display.max_columns', None)\n",
    "from pandas.tseries.offsets import MonthEnd # 月末\n",
    "Market_ret = pd.read_csv('Marketret_mon_stock2023.csv')\n",
    "Market_ret['month'] = pd.to_datetime(Market_ret['month'], format='%b %Y') + MonthEnd(0)\n",
    "Market_ret.set_index('month', inplace=True)\n",
    "Market_ret.sort_index(inplace=True)\n",
    "Market_ret            \n",
    "inflation = pd.read_csv('inflation.csv')\n",
    "inflation['month'] = pd.to_datetime(inflation['month'],format='%Y/%m/%d')\n",
    "inflation.set_index('month',inplace=True)\n",
    "inflation\n",
    "price_dividend = pd.read_csv('Price_dividend_mon2024.csv')\n",
    "price_dividend['month'] = pd.date_range(start='1990-12-31', end='2024-12-31', freq='ME')\n",
    "price_dividend.set_index('month', inplace=True)\n",
    "price_dividend.sort_index(inplace=True)\n",
    "price_dividend = price_dividend.drop(columns=['Unnamed: 0'])\n",
    "\n",
    "price_earning = pd.read_csv('Price_earnings_mon2024.csv')\n",
    "price_earning['month'] = pd.date_range(start='1991-01-31', end='2024-12-31', freq='ME')\n",
    "price_earning.set_index('month', inplace=True)\n",
    "price_earning.sort_index(inplace=True)\n",
    "\n",
    "price_bookvalue = pd.read_csv('Price_bookvalue_mon2024.csv')\n",
    "price_bookvalue['month'] = pd.date_range(start='1990-12-31', end='2024-12-31', freq='ME')\n",
    "price_bookvalue.set_index('month', inplace=True)\n",
    "price_bookvalue.sort_index(inplace=True)\n",
    "Market_ret_day = pd.read_csv('Marketret_day_stock2024 .csv')\n",
    "Market_ret_day[\"Day\"] = pd.to_datetime(Market_ret_day[\"Day\"],format='%Y/%m/%d')\n",
    "Market_ret_day.set_index('Day',inplace=True)\n",
    "Market_ret_day.sort_index(inplace=True)\n",
    "Market_ret_day\n",
    "# Calculate monthly realized variance (RV) by summing squared excess returns\n",
    "Market_variance = Market_ret_day.resample('ME').apply(lambda df: pd.Series({\n",
    "    'RV': (df['er']**2).sum(),\n",
    "    'RV1': (df['er']**2).sum() + 2*(df['er']*df['er1']).sum(),\n",
    "    'RV2': (df['er']**2).sum() + 2*(df['er']*df['er1']).sum() + 2*(df['er']*df['er2']).sum(),\n",
    "    'RV3': (df['er']**2).sum() + 2*(df['er']*df['er1']).sum() + 2*(df['er']*df['er2']).sum() + 2*(df['er']*df['er3']).sum()\n",
    "}))\n",
    "\n",
    "Market_variance['var'] = Market_ret_day.resample('ME')['er'].var()\n",
    "Market_variance.index.name = 'month'\n",
    "\n",
    "# if RV3 <0, set to RV2 if RV2 <0, set to RV1 if RV1 <0, set to RV\n",
    "Market_variance['RV1'] = Market_variance['RV1'].where(Market_variance['RV1']>=0, Market_variance['RV'])\n",
    "Market_variance['RV2'] = Market_variance['RV2'].where(Market_variance['RV2']>=0, Market_variance['RV1'])\n",
    "Market_variance['RV3'] = Market_variance['RV3'].where(Market_variance['RV3']>=0, Market_variance['RV2'])\n",
    "Market_variance\n",
    "\n",
    "# market_variance <- daily_data[,.(MV = sum(er^2)),by = 'month'] This is R code\n",
    "reg_data = pd.merge(Market_ret,inflation,on = 'month')\n",
    "reg_data = pd.merge(reg_data,Market_variance,on = 'month')\n",
    "reg_data = pd.merge(reg_data,price_dividend,on = 'month')\n",
    "reg_data = pd.merge(reg_data,price_earning,on = 'month')\n",
    "reg_data = pd.merge(reg_data,price_bookvalue,on = 'month')\n",
    "reg_data = reg_data[['MarketR','rfmonth','ret','cpi','RV','RV1','RV2','RV3','var','pd','pe','pb','marketret3','marketret6','marketret12']]\n",
    "# Output reg_data to reg_data.csv\n",
    "reg_data.to_csv('C:/Users/ma970/Desktop/python-and-data-analysis-master/python-and-data-analysis-master/datasets/reg_data.csv')\n",
    "# Output reg_data to reg_data.xlsx\n",
    "reg_data.to_excel('C:/Users/ma970/Desktop/python-and-data-analysis-master/python-and-data-analysis-master/datasets/reg_data.xlsx')\n",
    "\n",
    "reg_data\n",
    "# 对比RV3和RV var\n",
    "reg_data_plot = reg_data['2000-01':'2024-12'].copy()\n",
    "# Plot the RV3 and var into one graph\n",
    "fig, ax1 = plt.subplots(figsize=(10,5))\n",
    "\n",
    "ax1.plot(reg_data_plot['RV3'],color='blue',marker='o',linewidth=0.8,\n",
    "         markersize=4,\n",
    "         linestyle='-',label='Realized Variance (RV3)')\n",
    "ax1.set_ylabel('Realized Variance (RV3)',color='blue')\n",
    "ax1.tick_params(axis='y', labelcolor='blue')   \n",
    "ax1.set_title(\"China Stock Market Realized Variance (RV3) and Variance (var)\", fontsize=16)\n",
    "\n",
    "ax1.plot(reg_data_plot['var'],color='green',marker='o',linewidth=0.8,\n",
    "         markersize=4,\n",
    "         linestyle='-',label='Variance (var)')\n",
    "ax1.plot(reg_data_plot['RV'],color='red',marker='o',linewidth=0.8,\n",
    "         markersize=4,\n",
    "         linestyle='-',label='Realized Variance (RV)')\n",
    "\n",
    "ax1.legend(loc='upper left')\n",
    "plt.show();\n",
    "reg_data_plot = reg_data['2000-01':'2024-12'].copy()\n",
    "# Plot the China's stock market return and inflation into one graph\n",
    "fig, ax1 = plt.subplots(figsize=(10,5))\n",
    "# the linewidth and marker size are set to be very small\n",
    "ax1.plot(reg_data_plot['ret'],color='red',marker='o',linewidth=0.8,\n",
    "         markersize=4,\n",
    "         linestyle='-',label='China Stock Market Return')\n",
    "ax1.set_ylabel('China Stock Market Return',color='red')\n",
    "#ax1.set_xlabel('Month')\n",
    "\n",
    "# 设置x轴的日期显示格式\n",
    "data_format = mdates.DateFormatter('%Y')\n",
    "ax1.xaxis.set_major_formatter(data_format)\n",
    "ax1.xaxis.set_major_locator(mdates.YearLocator())\n",
    "\n",
    "# 转置x轴的日期显示格式\n",
    "plt.xticks(rotation = 90)\n",
    "\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(reg_data_plot['RV3'].shift(1),color='blue',marker='o',linewidth=0.8,\n",
    "         markersize=4,\n",
    "         linestyle='-',label='Realized Variance')\n",
    "\n",
    "ax2.set_ylabel('Realized Variance',color='blue')\n",
    "\n",
    "plt.title('China Stock Market Return and Realized Variance')\n",
    "\n",
    "# change the legend into one box\n",
    "lines, labels = ax1.get_legend_handles_labels()\n",
    "lines2, labels2 = ax2.get_legend_handles_labels()\n",
    "ax2.legend(lines + lines2, labels + labels2, loc='upper right')\n",
    "\n",
    "plt.show();\n",
    "reg_data['RV'].describe().round(5)\n",
    "reg_data['RV'].skew()\n",
    "reg_data['RV'].kurt()\n",
    "from statsmodels.tsa.stattools import adfuller as ADF\n",
    "\n",
    "# 对月收益率数据进行ADF检验\n",
    "adf_result = ADF(reg_data[reg_data.index >= '2000-01-31']['RV'])\n",
    "\n",
    "print('原始序列的ADF检验结果:')\n",
    "print(f'ADF Statistic: {adf_result[0]:.4f}')\n",
    "print(f'p-value: {adf_result[1]:.4f}')\n",
    "print('Critical Values:')\n",
    "for key, value in adf_result[4].items():\n",
    "    print(f'   {key}: {value:.4f}')\n",
    "\n",
    "if adf_result[1] <= 0.05:\n",
    "    print('结论: p-value小于0.05，拒绝原假设，序列是平稳的。')\n",
    "else:\n",
    "    print('结论: p-value大于0.05，未能拒绝原假设，序列是非平稳的。')\n",
    "reg_data['lRV'] = reg_data['RV'].shift(1)\n",
    "model_cpi = smf.ols('ret ~ lRV',\n",
    "                 data=reg_data['1995-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 6})\n",
    "print(model_cpi.summary())\n",
    "reg_data['lcpi'] = reg_data['cpi'].shift(2)\n",
    "model_cpi = smf.ols('ret ~ lcpi',\n",
    "                 data=reg_data['2000-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 6})\n",
    "print(model_cpi.summary())\n",
    "model_twovariables = smf.ols('ret ~ lRV + lcpi',\n",
    "                 data=reg_data['2000-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 6})\n",
    "print(model_twovariables.summary())\n",
    "reg_data['lpd'] = reg_data['pd'].shift(1)\n",
    "model3 = smf.ols('ret ~ lRV + lpd + lcpi',\n",
    "                 data=reg_data['2000-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 6})\n",
    "print(model3.summary())\n",
    "# Calculate monthly realized variance (RV) by summing squared excess returns\n",
    "Market_variance_Q = Market_ret_day.resample('QE').apply(lambda df: pd.Series({\n",
    "    'RV': (df['er']**2).sum(),\n",
    "    'RV1': (df['er']**2).sum() + 2*(df['er']*df['er1']).sum(),\n",
    "    'RV2': (df['er']**2).sum() + 2*(df['er']*df['er1']).sum() + 2*(df['er']*df['er2']).sum(),\n",
    "    'RV3': (df['er']**2).sum() + 2*(df['er']*df['er1']).sum() + 2*(df['er']*df['er2']).sum() + 2*(df['er']*df['er3']).sum()\n",
    "}))\n",
    "\n",
    "Market_variance_Q['var'] = Market_ret_day.resample('ME')['er'].var()\n",
    "Market_variance_Q.index.name = 'month'\n",
    "\n",
    "# if RV3 <0, set to RV2 if RV2 <0, set to RV1 if RV1 <0, set to RV\n",
    "Market_variance_Q['RV1'] = Market_variance_Q['RV1'].where(Market_variance_Q['RV1']>=0, Market_variance_Q['RV'])\n",
    "Market_variance_Q['RV2'] = Market_variance_Q['RV2'].where(Market_variance_Q['RV2']>=0, Market_variance_Q['RV1'])\n",
    "Market_variance_Q['RV3'] = Market_variance_Q['RV3'].where(Market_variance_Q['RV3']>=0, Market_variance_Q['RV2'])\n",
    "Market_variance_Q.index.name = 'Q'\n",
    "Market_variance_Q\n",
    "# Market_ret_day复合成季收益率数据\n",
    "Q_marketret = reg_data['1995-01':'2024-12'].resample('QE').apply(\n",
    "    lambda df: pd.Series({\n",
    "        'MarketR': np.prod(1 + df['MarketR']) - 1,\n",
    "        'rfqtr': np.prod(1 + df['rfmonth']) - 1,\n",
    "        'cpi': sum(df['cpi']),\n",
    "        'pd': df['pd'].iloc[-1],\n",
    "        'pe': df['pe'].iloc[-1],\n",
    "        'pb': df['pb'].iloc[-1]\n",
    "    })\n",
    ")\n",
    "Q_marketret['ret'] = Q_marketret['MarketR'] - Q_marketret['rfqtr']\n",
    "Q_marketret.index.name = 'Q'\n",
    "Q_marketret\n",
    "Qreg_data = pd.merge(Q_marketret,Market_variance_Q,on = 'Q')\n",
    "Qreg_data\n",
    "# Plot the China's stock market return and inflation into one graph\n",
    "fig, ax1 = plt.subplots(figsize=(10,5))\n",
    "# the linewidth and marker size are set to be very small\n",
    "ax1.plot(Qreg_data['ret'],color='red',marker='o',linewidth=0.8,\n",
    "         markersize=4,\n",
    "         linestyle='-',label='China Stock Market Return')\n",
    "ax1.set_ylabel('China Stock Market Return',color='red')\n",
    "#ax1.set_xlabel('Month')\n",
    "\n",
    "# 设置x轴的日期显示格式\n",
    "data_format = mdates.DateFormatter('%Y')\n",
    "ax1.xaxis.set_major_formatter(data_format)\n",
    "ax1.xaxis.set_major_locator(mdates.YearLocator())\n",
    "\n",
    "# 转置x轴的日期显示格式\n",
    "plt.xticks(rotation = 90)\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(Qreg_data['RV3'].shift(1),color='blue',marker='o',linewidth=0.8,\n",
    "         markersize=4,\n",
    "         linestyle='-',label='Realized Variance')\n",
    "\n",
    "ax2.set_ylabel('Realized Variance',color='blue')\n",
    "\n",
    "plt.title('China Stock Market Return and Realized Variance')\n",
    "\n",
    "# change the legend into one box\n",
    "lines, labels = ax1.get_legend_handles_labels()\n",
    "lines2, labels2 = ax2.get_legend_handles_labels()\n",
    "ax2.legend(lines + lines2, labels + labels2, loc='upper right')\n",
    "\n",
    "# save figure\n",
    "import os\n",
    "# 检查并创建images目录\n",
    "if not os.path.exists('images'):\n",
    "    os.makedirs('images')\n",
    "fig.savefig('images/China Stock Market Return and Realized Variance Quarter.png',dpi = 1000,bbox_inches='tight')\n",
    "\n",
    "plt.show();\n",
    "Qreg_data['lRV'] = Qreg_data['RV'].shift(1)\n",
    "Qreg_data['lRV1'] = Qreg_data['RV1'].shift(1)\n",
    "Qreg_data['lRV2'] = Qreg_data['RV2'].shift(1)\n",
    "Qreg_data['lRV3'] = Qreg_data['RV3'].shift(1)\n",
    "model_qrv = smf.ols('ret ~ lRV',\n",
    "                 data=Qreg_data['1995-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model_qrv1 = smf.ols('ret ~ lRV1',\n",
    "                 data=Qreg_data['1995-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model_qrv2 = smf.ols('ret ~ lRV2',\n",
    "                 data=Qreg_data['1995-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model_qrv3 = smf.ols('ret ~ lRV3',\n",
    "                 data=Qreg_data['1995-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "\n",
    "# print all model summaries together\n",
    "print(\"=\" * 80)\n",
    "print(\"Model with lRV:\")\n",
    "print(\"=\" * 80)\n",
    "print(model_qrv.summary())\n",
    "print(\"\\n\" + \"=\" * 80)\n",
    "print(\"Model with lRV1:\")\n",
    "print(\"=\" * 80)\n",
    "print(model_qrv1.summary())\n",
    "print(\"\\n\" + \"=\" * 80)\n",
    "print(\"Model with lRV2:\")\n",
    "print(\"=\" * 80)\n",
    "print(model_qrv2.summary())\n",
    "print(\"\\n\" + \"=\" * 80)\n",
    "print(\"Model with lRV3:\")\n",
    "print(\"=\" * 80)\n",
    "print(model_qrv3.summary())\n",
    "# print all model summaries together\n",
    "!pip install stargazer\n",
    "from stargazer.stargazer import Stargazer\n",
    "from IPython.core.display import HTML\n",
    "\n",
    "stargazer = Stargazer([model_qrv, model_qrv1, model_qrv2, model_qrv3])\n",
    "stargazer.title(\"回归结果对比分析\")\n",
    "stargazer.custom_columns(['模型一', '模型二', '模型三', '模型四'])\n",
    "stargazer.show_model_numbers(False)\n",
    "stargazer.show_confidence_intervals(False)\n",
    "\n",
    "# 生成HTML格式\n",
    "html_output = stargazer.render_html()\n",
    "HTML(html_output)\n",
    "# 使用 summary_col 进行多模型对比（推荐用于Jupyter）\n",
    "from statsmodels.iolib.summary2 import summary_col\n",
    "\n",
    "# 创建对比表\n",
    "results_table = summary_col(\n",
    "    [model_qrv, model_qrv1, model_qrv2, model_qrv3],\n",
    "    model_names=['RV', 'RV1', 'RV2', 'RV3'],\n",
    "    stars=True,\n",
    "    float_format='%.3f',\n",
    "    info_dict={\n",
    "        'N': lambda x: f\"{int(x.nobs)}\"\n",
    "    }\n",
    ")\n",
    "\n",
    "print(\"\\n 使用 statsmodels 的 summary_col（显示t值）\\n\")\n",
    "print(results_table)\n",
    "Qreg_data['lpd'] = Qreg_data['pd'].shift(1)\n",
    "Qreg_data['lpb'] = Qreg_data['pb'].shift(1)\n",
    "Qreg_data['lpe'] = Qreg_data['pe'].shift(1)\n",
    "model1 = smf.ols('ret ~ lRV + lpb',\n",
    "                    data=Qreg_data['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model2 = smf.ols('ret ~ lRV + lpe',\n",
    "                    data=Qreg_data['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model3 = smf.ols('ret ~ lRV + lpd',\n",
    "                    data=Qreg_data['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "\n",
    "# 创建对比表\n",
    "results_table = summary_col(\n",
    "    [model1, model2, model3],\n",
    "    model_names=['PB模型', 'PE模型', 'PD模型'],\n",
    "    stars=True,\n",
    "    float_format='%.3f',\n",
    "    info_dict={\n",
    "        'N': lambda x: f\"{int(x.nobs)}\"\n",
    "    }\n",
    ")\n",
    "print(results_table)\n",
    "model = smf.ols('ret ~ lRV3',\n",
    "                    data=Qreg_data['1995-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model1 = smf.ols('ret ~ lRV3 + lpb',\n",
    "                    data=Qreg_data['1995-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model2 = smf.ols('ret ~ lRV3 + lpe',\n",
    "                    data=Qreg_data['1995-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model3 = smf.ols('ret ~ lRV3 + lpd',\n",
    "                    data=Qreg_data['1995-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "\n",
    "# 创建对比表\n",
    "results_table = summary_col(\n",
    "    [model, model1, model2, model3],\n",
    "    model_names=['RV3', 'PB', 'PE', 'PD'],\n",
    "    stars=True,\n",
    "    float_format='%.3f',\n",
    "    info_dict={\n",
    "        'N': lambda x: f\"{int(x.nobs)}\"\n",
    "    }\n",
    ")\n",
    "print(results_table)\n",
    "Qreg_data['lcpi'] = Qreg_data['cpi'].shift(1)\n",
    "model = smf.ols('ret ~ lRV3',\n",
    "                    data=Qreg_data['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model1 = smf.ols('ret ~ lRV3 + lpb',\n",
    "                    data=Qreg_data['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model2 = smf.ols('ret ~ lRV3 + lpe',\n",
    "                    data=Qreg_data['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model3 = smf.ols('ret ~ lRV3 + lpd',\n",
    "                    data=Qreg_data['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "\n",
    "model4 = smf.ols('ret ~ lRV3 + lcpi',\n",
    "                    data=Qreg_data['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model5 = smf.ols('ret ~ lRV3 + lpb + lcpi',\n",
    "                    data=Qreg_data['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model6 = smf.ols('ret ~ lRV3 + lpe + lcpi',\n",
    "                    data=Qreg_data['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model7 = smf.ols('ret ~ lRV3 + lpd + lcpi',\n",
    "                    data=Qreg_data['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "\n",
    "# 创建对比表\n",
    "results_table = summary_col(\n",
    "    [model, model1, model2, model3, model4, model5, model6, model7],\n",
    "    model_names=['RV3', 'PB', 'PE', 'PD', 'CPI', 'PB+CPI', 'PE+CPI', 'PD+CPI'],\n",
    "    stars=True,\n",
    "    float_format='%.3f',\n",
    "    info_dict={\n",
    "        'N': lambda x: f\"{int(x.nobs)}\"\n",
    "    }\n",
    ")\n",
    "print(results_table)\n",
    "# GARCH-in-Mean 模型拟合（使用 ARCHInMean）\n",
    "!pip install arch\n",
    "from arch.univariate import ARCHInMean, GARCH, Normal\n",
    "\n",
    "# 准备数据：选择时间段并转换为百分比单位\n",
    "garch_data = Market_ret_day.loc['2000-01-01':'2024-12-31', ['er']].copy() * 100\n",
    "garch_data = garch_data.dropna()\n",
    "\n",
    "# 使用 ARCHInMean 构建 GARCH(1,1)-in-Mean 模型\n",
    "print(\"=\"*80)\n",
    "print(\"GARCH(1,1)-in-Mean 模型估计\")\n",
    "print(\"=\"*80)\n",
    "\n",
    "# 创建 GARCH-in-Mean 模型\n",
    "# ARCHInMean 是专门用于 GARCH-in-Mean 的均值模型\n",
    "model_garch_m = ARCHInMean(garch_data['er'], volatility=GARCH())\n",
    "\n",
    "res = model_garch_m.fit(disp='off')\n",
    "\n",
    "print(res.summary())\n",
    "# Visualization: Conditional Volatility and Returns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Extract conditional volatility and variance from fitted results\n",
    "conditional_volatility = res.conditional_volatility\n",
    "conditional_variance = conditional_volatility\n",
    "returns = garch_data['er']\n",
    "\n",
    "# Create three subplots\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 10), sharex=True)\n",
    "\n",
    "# First plot: Returns\n",
    "ax1.plot(returns.index, returns.values, color='blue', linewidth=0.5, alpha=0.7, label='Excess Returns')\n",
    "ax1.axhline(y=0, color='black', linestyle='--', linewidth=0.5)\n",
    "ax1.set_ylabel('Daily Excess Return (%)', fontsize=11)\n",
    "ax1.set_title('Chinese Stock Market: GARCH-in-Mean Model Analysis', fontsize=14, fontweight='bold')\n",
    "ax1.grid(True, alpha=0.3)\n",
    "ax1.legend(loc='upper left')\n",
    "\n",
    "# Second plot: Conditional Volatility (Standard Deviation)\n",
    "ax2.plot(conditional_volatility.index, conditional_volatility.values, \n",
    "         color='red', linewidth=0.8, label='Conditional Volatility (σ_t)')\n",
    "ax2.set_ylabel('Conditional Volatility (%)', fontsize=11)\n",
    "ax2.grid(True, alpha=0.3)\n",
    "ax2.legend(loc='upper left')\n",
    "\n",
    "# Third plot: Conditional Variance\n",
    "ax3.plot(conditional_variance.index, conditional_variance.values, \n",
    "         color='green', linewidth=0.8, label='Conditional Variance (σ²_t)')\n",
    "ax3.set_ylabel('Conditional Variance (%²)', fontsize=11)\n",
    "ax3.set_xlabel('Date', fontsize=11)\n",
    "ax3.grid(True, alpha=0.3)\n",
    "ax3.legend(loc='upper left')\n",
    "\n",
    "# Mark important events on all subplots\n",
    "for ax in [ax1, ax2, ax3]:\n",
    "    # 2008 Financial Crisis\n",
    "    ax.axvspan(pd.Timestamp('2008-01-01'), pd.Timestamp('2009-01-01'), \n",
    "               alpha=0.15, color='gray')\n",
    "    # 2015 Stock Market Crash\n",
    "    ax.axvspan(pd.Timestamp('2015-06-01'), pd.Timestamp('2015-09-01'), \n",
    "               alpha=0.15, color='orange')\n",
    "    # 2020 COVID-19 Pandemic\n",
    "    ax.axvspan(pd.Timestamp('2020-01-01'), pd.Timestamp('2020-04-01'), \n",
    "               alpha=0.15, color='purple')\n",
    "\n",
    "# Add legend for historical events on the first plot\n",
    "from matplotlib.patches import Patch\n",
    "legend_elements = [\n",
    "    Patch(facecolor='gray', alpha=0.3, label='2008 Financial Crisis'),\n",
    "    Patch(facecolor='orange', alpha=0.3, label='2015 Market Crash'),\n",
    "    Patch(facecolor='purple', alpha=0.3, label='2020 COVID-19')\n",
    "]\n",
    "ax1.legend(handles=legend_elements, loc='upper right', fontsize=9)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('images/GARCH_M_Complete_Analysis.png', dpi=300, bbox_inches='tight')\n",
    "plt.show();\n",
    "# GARCH-in-Mean Model with Monthly Data\n",
    "from arch.univariate import ARCHInMean, GARCH, Normal\n",
    "reg_data['lpb'] = reg_data['pb'].shift(1)\n",
    "# Prepare monthly data: select time period and convert to percentage\n",
    "reg_mon = reg_data.loc['2000-01':'2024-12', ['ret','lpb']].copy() \n",
    "reg_mon['ret'] = reg_mon['ret'] * 100\n",
    "reg_mon = reg_mon.dropna()\n",
    "\n",
    "# Build GARCH(1,1)-in-Mean model with monthly data\n",
    "print(\"=\"*80)\n",
    "print(\"GARCH(1,1)-in-Mean Model Estimation (Monthly Data)\")\n",
    "print(\"=\"*80)\n",
    "\n",
    "# Create GARCH-in-Mean model\n",
    "model_garch_m_mon = ARCHInMean(y=reg_mon['ret'],x=reg_mon['lpb'], volatility=GARCH())\n",
    "\n",
    "res_mon = model_garch_m_mon.fit(disp='off')\n",
    "\n",
    "print(res_mon.summary())\n",
    "# Visualization: Monthly GARCH-in-Mean Results\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Extract conditional volatility and variance from fitted results\n",
    "conditional_volatility_mon = res_mon.conditional_volatility\n",
    "conditional_variance_mon = conditional_volatility_mon**2\n",
    "returns_mon = reg_mon['ret']\n",
    "\n",
    "# Create three subplots\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 10), sharex=True)\n",
    "\n",
    "# First plot: Returns\n",
    "ax1.plot(returns_mon.index, returns_mon.values, color='blue', linewidth=0.8, \n",
    "         marker='o', markersize=3, alpha=0.7, label='Excess Returns')\n",
    "ax1.axhline(y=0, color='black', linestyle='--', linewidth=0.5)\n",
    "ax1.set_ylabel('Monthly Excess Return (%)', fontsize=11)\n",
    "ax1.set_title('Chinese Stock Market: GARCH-in-Mean Model (Monthly Data)', fontsize=14, fontweight='bold')\n",
    "ax1.grid(True, alpha=0.3)\n",
    "ax1.legend(loc='upper left')\n",
    "\n",
    "# Second plot: Conditional Volatility (Standard Deviation)\n",
    "ax2.plot(conditional_volatility_mon.index, conditional_volatility_mon.values, \n",
    "         color='red', linewidth=1.0, label='Conditional Volatility (σ_t)')\n",
    "ax2.set_ylabel('Conditional Volatility (%)', fontsize=11)\n",
    "ax2.grid(True, alpha=0.3)\n",
    "ax2.legend(loc='upper left')\n",
    "\n",
    "# Third plot: Conditional Variance\n",
    "ax3.plot(conditional_variance_mon.index, conditional_variance_mon.values, \n",
    "         color='green', linewidth=1.0, label='Conditional Variance (σ²_t)')\n",
    "ax3.set_ylabel('Conditional Variance (%²)', fontsize=11)\n",
    "ax3.set_xlabel('Date', fontsize=11)\n",
    "ax3.grid(True, alpha=0.3)\n",
    "ax3.legend(loc='upper left')\n",
    "\n",
    "# Mark important events on all subplots\n",
    "for ax in [ax1, ax2, ax3]:\n",
    "    # 2008 Financial Crisis\n",
    "    ax.axvspan(pd.Timestamp('2008-01-01'), pd.Timestamp('2009-01-01'), \n",
    "               alpha=0.15, color='gray')\n",
    "    # 2015 Stock Market Crash\n",
    "    ax.axvspan(pd.Timestamp('2015-06-01'), pd.Timestamp('2015-09-01'), \n",
    "               alpha=0.15, color='orange')\n",
    "    # 2020 COVID-19 Pandemic\n",
    "    ax.axvspan(pd.Timestamp('2020-01-01'), pd.Timestamp('2020-04-01'), \n",
    "               alpha=0.15, color='purple')\n",
    "\n",
    "# Add legend for historical events on the first plot\n",
    "from matplotlib.patches import Patch\n",
    "legend_elements = [\n",
    "    Patch(facecolor='gray', alpha=0.3, label='2008 Financial Crisis'),\n",
    "    Patch(facecolor='orange', alpha=0.3, label='2015 Market Crash'),\n",
    "    Patch(facecolor='purple', alpha=0.3, label='2020 COVID-19')\n",
    "]\n",
    "ax1.legend(handles=legend_elements, loc='upper right', fontsize=9)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('images/GARCH_M_Monthly_Analysis.png', dpi=300, bbox_inches='tight')\n",
    "plt.show();\n",
    "reg_data['lRV3'] = reg_data['RV3'].shift(1)\n",
    "model_rv = smf.ols('RV3~ lRV3',\n",
    "                 data=reg_data['2000-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 6})\n",
    "print(model_rv.summary())\n",
    "Turnover = pd.read_csv('Turnover_individual_mon2022.csv')\n",
    "Turnover['month'] = pd.to_datetime(Turnover['month'])\n",
    "Turnover.set_index('month', inplace=True)\n",
    "if len(Turnover) <= 1000:  # 控制生成日期的数量\n",
    "    Turnover['month'] = pd.date_range(start='1990-12-31', periods=len(Turnover), freq='M')\n",
    "    Turnover.set_index('month', inplace=True)\n",
    "else:\n",
    "    print(\"数据长度过大，无法生成有效日期\")\n",
    "Turnover.sort_index(inplace=True)\n",
    "Turnover\n",
    "fig, ax1 = plt.subplots(figsize=(10,5))\n",
    "ax1.plot(Turnover['to_v'],color='blue',marker='o',linewidth=0.8,\n",
    "         markersize=4,\n",
    "         linestyle='-',label='Market Turnover Ratio')\n",
    "ax1.set_ylabel('Market Turnover Ratio',color='blue')\n",
    "#ax1.set_xlabel('Month')\n",
    "data_format = mdates.DateFormatter('%Y')\n",
    "ax1.xaxis.set_major_formatter(data_format)\n",
    "ax1.xaxis.set_major_locator(mdates.YearLocator())\n",
    "plt.xticks(rotation = 90)\n",
    "\n",
    "plt.title('China Stock Market Turnover Ratio', fontsize=16)\n",
    "plt.show();\n",
    "# RV and Turnover correlation\n",
    "reg_data_turnover = pd.merge(reg_data,Turnover,on='month',how='left')\n",
    "test = reg_data_turnover['2000-01':'2024-12']\n",
    "test[['RV3','to_v']].corr()\n",
    "fig, ax1 = plt.subplots(figsize=(10,5))\n",
    "ax1.plot(test['RV3'],color='blue',marker='o',linewidth=0.8,\n",
    "         markersize=4,\n",
    "         linestyle='-',label='Realized Variance (RV3)')\n",
    "ax1.set_ylabel('Realized Variance (RV3)',color='blue')\n",
    "ax1.tick_params(axis='y', labelcolor='blue')   \n",
    "ax1.set_title(\"China Stock Market Realized Variance (RV3) and Turnover Ratio\", fontsize=16)\n",
    "\n",
    "ax1.set_ylabel('Realized Variance (RV3)',color='blue')\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(test['to_v'],color='green',marker='o',linewidth=0.8,\n",
    "         markersize=4,\n",
    "         linestyle='-',label='Market Turnover Ratio')\n",
    "ax2.set_ylabel('Market Turnover Ratio',color='green')\n",
    "ax2.tick_params(axis='y', labelcolor='green')\n",
    "\n",
    "lines, labels = ax1.get_legend_handles_labels()\n",
    "lines2, labels2 = ax2.get_legend_handles_labels()\n",
    "ax2.legend(lines + lines2, labels + labels2, loc='upper left')\n",
    "plt.show();\n",
    "reg_data_turnover['lto_v'] = reg_data_turnover['to_v'].shift(1)\n",
    "model_turnover = smf.ols('RV3 ~ lRV3 + lto_v',\n",
    "                 data=reg_data_turnover['2000-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 6})\n",
    "print(model_turnover.summary())\n",
    "reg_data_turnover['CV'] = model_turnover.fittedvalues\n",
    "reg_data_turnover['lpb'] = reg_data_turnover['pb'].shift(1)\n",
    "model_cv = smf.ols('ret ~ CV',\n",
    "                 data=reg_data_turnover['2000-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 6})\n",
    "model_cv2 = smf.ols('ret ~ CV + lpb',\n",
    "                    data=reg_data_turnover['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 6})\n",
    "print(model_cv.summary())\n",
    "print(model_cv2.summary())\n",
    "Q_turnover = Turnover.resample('QE').sum()\n",
    "Q_turnover.index.name = 'Q'\n",
    "Q_turnover\n",
    "Qreg_data_turnover = pd.merge(Qreg_data, Q_turnover, on='Q', how='left')\n",
    "Qreg_data_turnover['lto_v'] = Qreg_data_turnover['to_v'].shift(1)\n",
    "model_qturnover = smf.ols('RV3 ~ lRV3 + lto_v',\n",
    "                 data=Qreg_data_turnover['2000-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "print(model_qturnover.summary())\n",
    "Qreg_data_turnover['CV'] = model_qturnover.fittedvalues\n",
    "Qreg_data_turnover['lpb'] = Qreg_data_turnover['pb'].shift(1)\n",
    "model_qcv = smf.ols('ret ~ CV',\n",
    "                 data=Qreg_data_turnover['2000-01':'2024-12']).fit(\n",
    "                     cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model_qcv2 = smf.ols('ret ~ CV + lpb',\n",
    "                    data=Qreg_data_turnover['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "model_qcv3 = smf.ols('ret ~ CV + lpb + lcpi',\n",
    "                    data=Qreg_data_turnover['2000-01':'2024-12']).fit(\n",
    "                        cov_type='HAC', cov_kwds={'maxlags': 2})\n",
    "print(model_qcv.summary())\n",
    "print(model_qcv2.summary())\n",
    "print(model_qcv3.summary()) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cae8a4b8",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
